This code implements all statistical analyses and creates figures for “Longitudinal development of brain iron is linked to cognition in youth” (Larsen, et al., 2019).
Relevant functions:
GAMM modeling
# This function:
# 1. executes the GAMM model,
# 2. sends output to the parametric bootstrap (if requested),
# 3. prints a regression table, and
# 4. sends the model to the visualizer for plotting.
gamm_model <- function(df, model_formula, this_label, smooth_var, int_var = NULL,group_var, pbootstrap = F, longPlot = F, model_test = T){
cat(sprintf("\n\n### Results for %s\n",this_label))
model_formula <- as.formula(model_formula)
if(!"exclusions" %in% colnames(df)) {
df$exclusions <- FALSE; #there is no exclusions column so make one that excludes none
}
g1<-gamm(model_formula,
data=df,
random = list(bblid =~ 1),
subset = exclusions == F)
if (model_test == T){
if (pbootstrap == T) {
#Send to bootstrap function
g1$pb<-pboot(g1)
#Print a table that shows the bootstrap outcome
print(g1$pb %>%
summary() %>%
.$test %>%
as.data.frame() %>%
kable(caption = sprintf("Parametric Bootstrap test for %s",this_label)) %>%
kable_styling(full_width = F, position = "left",bootstrap_options = c("striped"))
)
if (isTRUE(all.equal(g1$pb$bestmod,model_formula))) {
cat("The initial (more complicated) model is best")
g <- g1
# Refit the model without fixed df for the purpose of plotting
plot_formula <- as.formula(gsub(", fx = T","",deparse(model_formula)))
cat(deparse(plot_formula))
plotg<-g
plotg <- gamm(plot_formula,
data = df,
random = list(bblid =~1),
subset = exclusions == F)
} else {
cat("The simpler model is best")
cat(" refitting ")
g <-gamm(as.formula(g1$pb$bestmod),
data=df,
random = list(bblid =~ 1),
subset = exclusions == F)
#Again, this part not implemented (model not refit with unfixed df)
plot_formula <- as.formula(gsub("ti\\(",'te\\(',deparse(g$gam$formula)) %>% gsub(", fx = T", "", .))
plotg<-g
plotg <-gamm(plot_formula,
data=df,
random = list(bblid =~ 1),
subset = exclusions == F)
}
} else {
if (!is.null(int_var)) {
# We are not bootstrapping, but there is an interaction variable
s<-summary(g1$gam)
if (s$s.table[grep(x=rownames(s$s.table),pattern = int_var),"p-value"] <.05/4) {
#Checked if interaction is sig, if so keep in the model
g <- g1
plot_formula <- as.formula(gsub(", fx = T", "", deparse(model_formula)))
plotg <- g
} else {
#Interaction is not sig, remove from the model
cat("The simpler model is best")
thisResp <- as.character(g1$gam$terms[[2]])
theseVars <- attr(terms(model_formula),"term.labels")
new_formula <- reformulate(theseVars[0:(length(theseVars)-1)],response = thisResp)
g <-gamm(as.formula(new_formula),
data=df,
random = list(bblid =~ 1),
subset = exclusions == F)
plot_formula <- as.formula(gsub("ti\\(",'te\\(',deparse(g$gam$formula)) %>% gsub(", fx = T", "", .))
plotg<-gamm(as.formula(plot_formula),
data=df,
random = list(bblid =~ 1),
subset = exclusions == F)
}
} else {
#There is no interaction term, just plot.
g <- g1
plot_formula <- as.formula(gsub(", fx = T", "", deparse(model_formula)))
plotg <-gamm(as.formula(plot_formula),
data=df,
random = list(bblid =~ 1),
subset = exclusions == F)
}
}
} else {
g <- g1
plotg<-g
}
g$gam$data<-df %>%
filter(exclusions == F)
#Display model results:
s_tidytable<- tidy(g$gam)
p_tidytable <- tidy(g$gam,parametric = T)
snames = names(s_tidytable)
pnames = names(p_tidytable)
names(s_tidytable)=pnames
thisBIC <- BIC(g$lme)
numObs <- g$lme$dims$N
g$BIC <- thisBIC
print(concurvity(g$gam)%>%kable(caption = "convurvity")%>%kable_styling(full_width = F,bootstrap_options = "striped",position = "left"))
stattable <- rbind(p_tidytable,snames,s_tidytable) %>%
kable(caption = sprintf("Regression table from gamm in %s, BIC = %1.2f, obs = %d",this_label,thisBIC,numObs)) %>%
kable_styling(full_width = F, position = "left")
print(stattable)
write.csv(x = rbind(p_tidytable,snames,s_tidytable),file = sprintf("GAMM_table_%s.csv",this_label),row.names = F)
cat(sprintf("Regression table from gamm in %s, BIC = %1.2f, obs = %d\n",this_label,thisBIC,numObs),
file = sprintf("GAMM_table_%s.txt",this_label))
cat(sprintf("Bootstrap p value %1.5f",g1$pb$test["PBtest","p.value"]),
file = sprintf("GAMM_table_%s.txt",this_label),
append = T)
#Send final model to visualizer:
if (longPlot == T) {
g$pl <- longitudinal_plot(g,plabels = this_label)
} else{
if (s_tidytable$p.value[nrow(s_tidytable)]<1.05) {
g$pl <- visualize_models(plotg,plabels = this_label, smooth_var = smooth_var, int_var = int_var, group_var = group_var)
}
}
#Return result object
result <- g
return(result)
}
Bootstrap procedure
## Parametric bootstrap of likelihood ratio test for nested models
pboot <- function(modelobj){
numsims <- 1000
df <- modelobj$gam$model
thisResp <- as.character(modelobj$gam$terms[[2]])
f1 <- modelobj$gam$formula
theseVars <- attr(terms(f1),"term.labels")
f2 <- reformulate(theseVars[0:(length(theseVars)-1)],response = thisResp)
g1 <- gam(f1,data = df)
g2 <- gam(f2,data = df)
mat1 <- model.matrix(g1)
mat2 <- model.matrix(g2)
bblid<-df$bblid
y <- df[,thisResp]
m1 <- lmer(y ~ -1 + mat1 + (1|bblid))
m2 <- lmer(y ~ -1 + mat2 + (1|bblid))
refdist <- PBrefdist(m1, m2, nsim=numsims)#, cl=cl)
pb <- PBmodcomp(m1, m2, ref = refdist)
int_pval <- pb$test["PBtest","p.value"]
if (int_pval < .05/4) {
pb$bestmod <- f1
} else {
pb$bestmod <- f2
}
return(pb)
}
Visualize models
Load the data
# Load the datafile
dataFile <- "brain_iron_data.Rda"
dataTable <- readRDS(file=dataFile)
## Create the final sample
### Remove
### 1. Projects with incompatible data
dataTable <- dataTable %>%
filter(bblid != 110689)%>%
filter(!(ProjectName %in% c("DAY2_808799",
"FNDM1_810211",
"FNDM2_810211",
"NEFF_V2",
"NEFF_PILOT"
)
)
) #Not using data from these projects or participants with no scans.
### 2. Sequences with odd parameters
dataTable <- dataTable %>%
filter(sequence != 'B0map_v4_matchedFOV' & sequence != "B0map") %>% # Different voxel sizes
filter(sequence != 'B0map_onesizefitsall_v3_T2S') %>% # small number of scans with odd params
filter(!is.na(GOOD)&!is.na(Putamen)) # No r2* data available or no QC possible
#Print some info
summaryTable<-dataTable%>%group_by(bblid,sex)%>%summarize(n=n())
cat(sprintf('Initial sample includes:\n %d individuals aged %1.2f - %1.2f (M = %1.3f, SD = %.3f)\n (M/F = %d/%d)\n%d observations',
length(unique(dataTable$bblid)),
min(dataTable$age,na.rm=T),
max(dataTable$age[dataTable$visit==1],na.rm=T),
mean(dataTable$age[dataTable$visit==1],na.rm=T),
sd(dataTable$age[dataTable$visit==1],na.rm=T),
sum(summaryTable$sex=='male'),
sum(summaryTable$sex=='female'),
length(dataTable$bblid)))
## Initial sample includes:
## 1543 individuals aged 8.17 - 26.92 (M = 15.237, SD = 3.738)
## (M/F = 728/815)
## 2321 observations
### 3. Subjects that don't meet inclusion
cat(sprintf('\nHealth exclusions: %d subs, %d scans',
length(unique(dataTable$bblid[dataTable$ltnExcludev2==1])),
sum(dataTable$ltnExcludev2==1,na.rm = T)))
##
## Health exclusions: 309 subs, 419 scans
dataTable <- dataTable %>%
filter(ltnExcludev2 ==0 | is.na(ltnExcludev2)) # exclude unhealthy/medicated subs
### 4. Scans that don't pass QA
cat(sprintf('\nQC exclusions: %d scans fail (%d pass)',
sum(dataTable$GOOD==0,na.rm = T),
sum(dataTable$GOOD %in% c("1","2"),na.rm = T)))
##
## QC exclusions: 641 scans fail (1236 pass)
dataTable <- dataTable %>%
filter(GOOD %in% c("1","2")) # keep only subjects that pass QA
## Count visits
dataTable <- dataTable %>%
group_by(bblid) %>%
mutate(visitnum = min_rank(ScanAgeMonths)) %>%
ungroup()
dataTable$visitnum <- ordered(dataTable$visitnum)
summaryTable<-dataTable%>%group_by(bblid,sex)%>%summarize(n=n())
cat(sprintf('\nFinal sample includes:\n %d individuals aged %1.2f - %1.2f (M = %1.3f, SD = %.3f)\n (M/F = %d/%d)\n %d observations',
length(unique(summaryTable$bblid)),
min(dataTable$age,na.rm=T),
max(dataTable$age[dataTable$visit==1],na.rm=T),
mean(dataTable$age[dataTable$visit==1],na.rm=T),
sd(dataTable$age[dataTable$visit==1],na.rm=T),
sum(summaryTable$sex=='male'),
sum(summaryTable$sex=='female'),
length(dataTable$bblid)))
##
## Final sample includes:
## 922 individuals aged 8.17 - 26.92 (M = 15.146, SD = 3.724)
## (M/F = 436/486)
## 1236 observations
# Create sample description figure (Fig 1)
sortedBBLID<-dataTable %>% group_by(bblid) %>% summarise(m=min(age)) %>% arrange(m)
sortedBBLID$row<-as.numeric(rownames(sortedBBLID))
newDataTable<-dataTable%>%left_join(sortedBBLID%>%select(bblid,row),by="bblid")
fig1<-ggplot(data = newDataTable,aes(x=reorder(row,age,FUN = min),y=age,color=oSex)) +
geom_point(size=.5,alpha=.5) + geom_line(alpha=.5) +
scale_x_discrete(breaks = c(1,500,1000,length(sortedBBLID$bblid))) +
coord_flip(clip = "off") +
scale_color_manual(values = c("blue","red"),labels=c("Male","Female"))+
labs(y = "Age (years)",x="Participant",color="") +theme(legend.position = c(.8,.2))
print(fig1)

ggsave("Fig1.svg",device = "svg",plot = fig1,width = 85,height = 85,units = "mm",bg = "transparent")
## make data frame for participants that have cognition data.
behTable <- dataTable %>%
filter(!is.na(NAR_Overall_Efficiency)) %>%
filter(scan2cnbmonths<=6)
summaryTable<-behTable%>%group_by(bblid,sex)%>%summarize(n=n())
cat(sprintf('\nBehavior sample includes:\n %d individuals aged %1.2f - %1.2f (M = %1.3f, SD = %.3f)\n (M/F = %d/%d)\n %d observations',
length(unique(summaryTable$bblid)),
min(behTable$age,na.rm=T),
max(behTable$age[behTable$visit==1],na.rm=T),
mean(behTable$age[behTable$visit==1],na.rm=T),
sd(behTable$age[behTable$visit==1],na.rm=T),
sum(summaryTable$sex=='male'),
sum(summaryTable$sex=='female'),
length(behTable$bblid)))
##
## Behavior sample includes:
## 818 individuals aged 8.17 - 26.92 (M = 14.840, SD = 3.574)
## (M/F = 389/429)
## 1067 observations
Combat harmonization
## Combat harmonization
# First harmonize for age models
comtable <- dataTable %>%
select(Caudate,Putamen,Accumbens_Area,Pallidum,sequence,age,oSex,bblid,scanid,visitnum) %>%
drop_na()
batch <- as.character(comtable$sequence)
ctab <- t(data.matrix(comtable%>%select(Caudate,Putamen,Accumbens_Area,Pallidum)))
g<-gamm(Putamen ~ visitnum + s(age, k = 4, fx = T) + oSex + s(age, by = oSex, k = 4, fx = T) ,data=comtable,random=list(bblid=~1))
mod<- model.matrix(g$gam)
combatdata <- combat(ctab,batch,mod=mod, eb = F)
## [combat] Performing ComBat without empirical Bayes (L/S model)
## [combat] Found 4 batches
## [combat] Adjusting for 9 covariate(s) or covariate level(s)
## [combat] Standardizing Data across features
## [combat] Fitting L/S model
## [combat] Adjusting the Data
harmonized_data<-data.frame(t(combatdata$dat.combat))
colnames(harmonized_data)<- paste(colnames(harmonized_data),"h",sep = "_")
harmonized_data$bblid=comtable$bblid
harmonized_data$scanid=comtable$scanid
dataTable<- dataTable %>%
left_join(harmonized_data,by = c("bblid","scanid"))
# Now for Behavior models
comtable <- behTable %>%
select(Caudate,Putamen,Accumbens_Area,Pallidum,sequence,age,NAR_Overall_Efficiency,oSex,bblid,scanid) %>%
drop_na()
batch <- as.character(comtable$sequence)
ctab <- t(data.matrix(comtable%>%select(Caudate,Putamen,Accumbens_Area,Pallidum)))
g<-gamm(Putamen ~ s(age, k = 4, fx = T) + s(age, by = NAR_Overall_Efficiency, k = 4, fx = T) + oSex,data=comtable,random=list(bblid=~1))
mod<- model.matrix(g$gam)
combatdata <- combat(ctab,batch,mod=mod, eb = F)
## [combat] Performing ComBat without empirical Bayes (L/S model)
## [combat] Found 3 batches
## [combat] Adjusting for 8 covariate(s) or covariate level(s)
## [combat] Standardizing Data across features
## [combat] Fitting L/S model
## [combat] Adjusting the Data
harmonized_data<-data.frame(t(combatdata$dat.combat))
colnames(harmonized_data)<- paste(colnames(harmonized_data),"h",sep = "_")
harmonized_data$bblid=comtable$bblid
harmonized_data$scanid=comtable$scanid
behTable<- behTable %>%
left_join(harmonized_data,by = c("bblid","scanid"))
longTable <- dataTable %>%
gather(key = "ROI",value="R2s",Accumbens_Area_h,Putamen_h,Caudate_h,Pallidum_h)
rois <- c("Caudate_h","Putamen_h","Accumbens_Area_h","Pallidum_h")
ROI comparison
Analyses for figure 3
source("/Users/larsenb/Documents/Tools/RainCloudPlots-master/tutorial_R/R_rainclouds.R") #requires RainCloud package: https://github.com/RainCloudPlots/RainCloudPlots/blob/master/tutorial_R/R_rainclouds.R
longTable$ROI<-factor(longTable$ROI,levels = c("Caudate_h","Putamen_h","Accumbens_Area_h","Pallidum_h"),labels = c("Caudate","Putamen","Nucleus Accumbens","Pallidum"))
longTable$oldROI <- c("Caudate","Putamen","Accumbens_Area","Pallidum")
l<-lmerTest::lmer(R2s ~ sex+age+ROI+(1|bblid),data = longTable)
m<-psycho::get_means(fit = l,formula="ROI")
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, set emm_options(pbkrtest.limit = 4944) or larger,
## but be warned that this may result in large computation time and memory use.
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, set emm_options(lmerTest.limit = 4944) or larger,
## but be warned that this may result in large computation time and memory use.
emmeans_table<- kable(m) %>%
kable_styling(bootstrap_options = "striped",position = "left")
print(emmeans_table)
|
ROI
|
Mean
|
SE
|
df
|
CI_lower
|
CI_higher
|
|
Caudate
|
15.17269
|
0.0552106
|
Inf
|
15.06448
|
15.28090
|
|
Putamen
|
16.95906
|
0.0552106
|
Inf
|
16.85085
|
17.06727
|
|
Nucleus Accumbens
|
17.58602
|
0.0552106
|
Inf
|
17.47781
|
17.69423
|
|
Pallidum
|
22.33203
|
0.0552106
|
Inf
|
22.22382
|
22.44024
|
c <- psycho::get_contrasts(l,"ROI",adjust="bonf")
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, set emm_options(pbkrtest.limit = 4944) or larger,
## but be warned that this may result in large computation time and memory use.
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, set emm_options(lmerTest.limit = 4944) or larger,
## but be warned that this may result in large computation time and memory use.
contrast_table<- kable(c) %>%
kable_styling(bootstrap_options = "striped",position = "left")
print(contrast_table)
|
Contrast
|
Difference
|
SE
|
df
|
z
|
p
|
CI_lower
|
CI_higher
|
|
Caudate - Putamen
|
-1.7863742
|
0.0676308
|
Inf
|
-26.413634
|
0
|
-1.9648016
|
-1.607947
|
|
Caudate - Nucleus Accumbens
|
-2.4133326
|
0.0676308
|
Inf
|
-35.683948
|
0
|
-2.5917600
|
-2.234905
|
|
Caudate - Pallidum
|
-7.1593439
|
0.0676308
|
Inf
|
-105.859280
|
0
|
-7.3377713
|
-6.980917
|
|
Putamen - Nucleus Accumbens
|
-0.6269584
|
0.0676308
|
Inf
|
-9.270314
|
0
|
-0.8053857
|
-0.448531
|
|
Putamen - Pallidum
|
-5.3729697
|
0.0676308
|
Inf
|
-79.445646
|
0
|
-5.5513970
|
-5.194542
|
|
Nucleus Accumbens - Pallidum
|
-4.7460113
|
0.0676308
|
Inf
|
-70.175332
|
0
|
-4.9244386
|
-4.567584
|
write.csv(c,'EMMeans_contrast.csv',row.names = F)
fig3<-ggplot(data=longTable,aes(x = reorder(ROI,-R2s,FUN = mean,na.rm=T),y = R2s,fill=ROI,color = ROI)) +
geom_flat_violin(show.legend = F,position=position_nudge(x=.2),alpha = .4) +
labs(x=NULL,y = "R2* (1/sec)",font_size=font_size) +
geom_point(size=.25,alpha=.4,position = position_jitter(width = .1),show.legend = F) +
geom_boxplot(color = "black",fill=NA,outlier.shape = NA,width=.1,show.legend = F) +
scale_color_brewer(palette = "Set1") + scale_fill_brewer(palette = "Set1") +
scale_x_discrete(labels=c("Globus\nPallidus","Nucleus\nAccumbens","Putamen","Caudate"))
fig3

ggsave(plot=fig3,"Fig3.svg",device = "svg",width = 85,height = 85,units = "mm",bg = "transparent")
Developmental effects
Analyses for figure 4.
Fit nonlinear age effect and check for age*sex interaction (factor smooth interaction: s(age, by = oSex)) oSex is Sex specified as an ordered factor to appropriately test for the interaction (female against male reference)
model_formula <- "R2s ~ visitnum +oSex + oRace_bwo +s(age, k = 4, fx = T) + s(age, by = oSex, k = 4, fx = T) "
## This will perform model testing. If the age * sex interaction is significant, it is retained, otherwise it will return an age + sex model.
models <- longTable %>%
group_by(ROI)%>%
nest()%>%
mutate(results=purrr::pmap(list(data,model_formula,this_label = ROI, pbootstrap = T, longPlot = T),.f=gamm_model))
Results for Nucleus Accumbens
Parametric Bootstrap test for Nucleus Accumbens
|
|
stat
|
df
|
ddf
|
p.value
|
|
PBtest
|
1.9757334
|
NA
|
NA
|
0.5664336
|
|
Gamma
|
1.9757334
|
NA
|
NA
|
0.5587682
|
|
Bartlett
|
2.0591528
|
3
|
NA
|
0.5602205
|
|
F
|
0.6585778
|
3
|
3.064699
|
0.6292186
|
|
LRT
|
1.9757334
|
3
|
NA
|
0.5774586
|
The simpler model is best refitting
convurvity
|
|
para
|
s(age)
|
|
worst
|
0.7753364
|
0.1211547
|
|
observed
|
0.7753364
|
0.0825105
|
|
estimate
|
0.7753364
|
0.1032679
|
Regression table from gamm in Nucleus Accumbens, BIC = 5844.61, obs = 1235
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
17.4861600325314
|
0.149943218021533
|
116.618545761905
|
0
|
|
visitnum.L
|
-0.0239854377469602
|
0.231762288339197
|
-0.103491546958909
|
0.917589783885288
|
|
visitnum.Q
|
0.0510209690673487
|
0.15345724144422
|
0.33247677716072
|
0.739586109756318
|
|
oSex.L
|
0.288357604672805
|
0.118627761102667
|
2.43077675910316
|
0.0152091226580745
|
|
oRace_bwo.L
|
-0.412873371047967
|
0.197447871953272
|
-2.09104999189698
|
0.0367290336123828
|
|
oRace_bwo.Q
|
-0.00173326477297276
|
0.153024269604343
|
-0.011326731226715
|
0.990964612066546
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
3
|
3
|
41.5735818338385
|
9.49021159068033e-26
|
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Sig change: 8.17 - 17.12
Results for Putamen
Parametric Bootstrap test for Putamen
|
|
stat
|
df
|
ddf
|
p.value
|
|
PBtest
|
8.96187
|
NA
|
NA
|
0.0299700
|
|
Gamma
|
8.96187
|
NA
|
NA
|
0.0305345
|
|
Bartlett
|
8.96842
|
3
|
NA
|
0.0297137
|
|
F
|
2.98729
|
3
|
3.001097
|
0.1963283
|
|
LRT
|
8.96187
|
3
|
NA
|
0.0298022
|
The simpler model is best refitting
convurvity
|
|
para
|
s(age)
|
|
worst
|
0.7753364
|
0.1211547
|
|
observed
|
0.7753364
|
0.1106914
|
|
estimate
|
0.7753364
|
0.1032679
|
Regression table from gamm in Putamen, BIC = 3944.17, obs = 1235
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
17.0837970708726
|
0.0705330905798225
|
242.209676769216
|
0
|
|
visitnum.L
|
0.243981166893351
|
0.11501228096822
|
2.12134882326843
|
0.0340926608837701
|
|
visitnum.Q
|
0.0247822592470953
|
0.0775572597245464
|
0.319535003365415
|
0.749375291685797
|
|
oSex.L
|
-0.0293967888188757
|
0.0529130616975227
|
-0.555567715716817
|
0.578607767500187
|
|
oRace_bwo.L
|
-0.0980967255048573
|
0.0881472318610424
|
-1.11287358018797
|
0.265980881783258
|
|
oRace_bwo.Q
|
0.155103316728515
|
0.068270910066593
|
2.27188002294424
|
0.023266395689981
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999999
|
2.99999999999999
|
86.1106136043234
|
3.18238513963805e-52
|
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Sig change: 10.05 - 26.92
Results for Caudate
Parametric Bootstrap test for Caudate
|
|
stat
|
df
|
ddf
|
p.value
|
|
PBtest
|
0.8792533
|
NA
|
NA
|
0.8231768
|
|
Gamma
|
0.8792533
|
NA
|
NA
|
0.8341081
|
|
Bartlett
|
0.8676429
|
3
|
NA
|
0.8332285
|
|
F
|
0.2930844
|
3
|
2.980323
|
0.8298361
|
|
LRT
|
0.8792533
|
3
|
NA
|
0.8304314
|
The simpler model is best refitting
convurvity
|
|
para
|
s(age)
|
|
worst
|
0.7753364
|
0.1211547
|
|
observed
|
0.7753364
|
0.0674739
|
|
estimate
|
0.7753364
|
0.1032679
|
Regression table from gamm in Caudate, BIC = 4271.13, obs = 1235
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
15.269762421595
|
0.0776776347370658
|
196.578622318796
|
0
|
|
visitnum.L
|
0.176402467706872
|
0.110934580845815
|
1.59014859353954
|
0.112059127248107
|
|
visitnum.Q
|
-0.00577971594470104
|
0.0716797635000482
|
-0.0806324639268259
|
0.9357474078279
|
|
oSex.L
|
0.146585890971771
|
0.0658920384329252
|
2.22463736830646
|
0.0262870751028674
|
|
oRace_bwo.L
|
-0.0586540660610204
|
0.10959347172685
|
-0.535196715067203
|
0.592610896545729
|
|
oRace_bwo.Q
|
0.230240038824217
|
0.0849748445558743
|
2.70950820831246
|
0.00683221278725056
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
3
|
3
|
9.95646441276744
|
1.73004055047032e-06
|
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Sig change: 8.17 - 15.89
Results for Pallidum
Parametric Bootstrap test for Pallidum
|
|
stat
|
df
|
ddf
|
p.value
|
|
PBtest
|
13.395682
|
NA
|
NA
|
0.0019980
|
|
Gamma
|
13.395682
|
NA
|
NA
|
0.0021317
|
|
Bartlett
|
13.224946
|
3
|
NA
|
0.0041746
|
|
F
|
4.465227
|
3
|
2.981003
|
0.1261747
|
|
LRT
|
13.395682
|
3
|
NA
|
0.0038546
|
The initial (more complicated) model is bestR2s ~ visitnum + oSex + oRace_bwo + s(age, k = 4) + s(age, by = oSex, k = 4)
convurvity
|
|
para
|
s(age)
|
s(age):oSexfemale
|
|
worst
|
0.7767801
|
0.6767412
|
0.6713269
|
|
observed
|
0.7767801
|
0.5599106
|
0.5935324
|
|
estimate
|
0.7767801
|
0.5484832
|
0.5320300
|
Regression table from gamm in Pallidum, BIC = 4753.09, obs = 1235
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
22.4420835917782
|
0.096170646920529
|
233.356895377061
|
0
|
|
visitnum.L
|
0.16947438464697
|
0.149994146382322
|
1.12987332328953
|
0.258751152755715
|
|
visitnum.Q
|
-0.041772753271621
|
0.0997112395016137
|
-0.418937258030425
|
0.6753355508665
|
|
oSex.L
|
-0.0883740132128409
|
0.0752483876130313
|
-1.17443065580765
|
0.240451197393393
|
|
oRace_bwo.L
|
-0.00948850745249523
|
0.12515802852913
|
-0.0758122156764943
|
0.939580896461433
|
|
oRace_bwo.Q
|
0.0685702907730544
|
0.0971245985102329
|
0.706003338235986
|
0.480320522722184
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999999
|
2.99999999999999
|
84.6072181107431
|
2.50042565031731e-51
|
|
s(age):oSexfemale
|
3.00000000000001
|
3
|
4.46988263857325
|
0.00395226040160837
|
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Sig change: 8.17 - 25.92
Sig change: 8.17 - 22.21 
Modeling of continuous*continuous interactions tests two possible models.
Analyses for figures 5 and 6.
We compare an additive (main effects + interaction) bivariate smooth model and a varying coefficient model (details below). The best model is chosen via the smallest BIC. After the best interaction model is selected, the significance of the interaction is tested with a parametric bootstrap likelihood ratio test. This test compares the model with the interaction term against a simpler nested model with main effects only. If the interaction model is significantly better, we keep that model. If not, the final model is the simpler model with no interaction.
2D smooth with additional main effects and tensor product smooth, ti: ti(age) + ti(Cognition) + ti(age,Cognition)
From documentation: This model specifies a main effects + interaction structure such as: y ~ ti(x) + ti(z) + ti(x,z)
ti is the proper way of specifying an interaction term in the context of included main effect terms:
“This functional ANOVA decomposition is supported by ti terms, which produce tensor product interactions from which the main effects have been excluded, under the assumption that they will be included separately. For example the ~ ti(x) + ti(z) + ti(x,z) would produce the above main effects + interaction structure. This is much better than attempting the same thing with s or te terms representing the interactions (although mgcv does not forbid it). Technically ti terms are very simple: they simply construct tensor product bases from marginal smooths to which identifiability constraints (usually sum-to-zero) have already been applied: correct nesting is then automatic (as with all interactions in a GLM framework). See Wood (2017, section 5.6.3).”
Varying coefficient model (using by =)
This will make the fit linear (rather than non-linear smooth) in the by variable (which is cognition in this case). From documentation: “When using by with a numberic covariate,”the by argument ensures that the smooth function gets multiplied by covariate z"
g<-gamm(NAR_Overall_Efficiency ~ s(age),data = behTable,random = list(bblid=~1))
r <- residuals(g$gam,type = "pearson")
behTable$Overall_Efficiency_resid <- r
g<-gamm(Putamen_h ~ s(age,k=4,fx=T),data = behTable,random = list(bblid=~1))
r <- residuals(g$gam,type = "pearson")
behTable$Putamen_h_resid <- r
# Loop over ROIs
for (r in "Putamen_h") {
cat(sprintf("\n## Results for %s\n",r))
# Set the cognitive variables to test
if (r == "Putamen_h") {
# Overall efficiency is significant for Putamen, so check the sub-domains of efficiency as well.
cog_vars <- c("NAR_Overall_Efficiency",
"NAR_F1_Social_Cognition_Efficiency",
"NAR_F2_Complex_Reasoning_Efficiency",
"NAR_F3_Executive_Efficiency",
"NAR_F4_Memory_Efficiency")
} else{
cog_vars <- c("NAR_Overall_Efficiency")
}
# Loop over the cognitive variables
for (cv in cog_vars) {
cat(sprintf("\n### %s\n",cv))
thisTable <- behTable
# First look at a main effects only model
# Fitting smooths for age and cognition
add_formula <- sprintf("%s ~ timepoint +oSex + s(age, k=4, fx = T) + s(%s, k=4, fx = T)",r,cv,cv)
gamm_model(thisTable,
add_formula,
this_label = sprintf("%s %s M.E.",r,cv),
smooth_var = cv,
group_var = "bblid",
pbootstrap = T,
model_test = F)
# Compare the two interaction models
cat('\n### Comparing interaction models...\n')
# Bivariate interaction
bv_formula <- sprintf(
"%s ~ timepoint +oSex + ti(age, k=4, fx = T) + ti(%s, k=4, fx = T) + ti(age,%s, k=4, fx = T)",
r,cv,cv,cv)
# Linear varying coefficient interaction
vc_formula <- sprintf(
"%s ~ timepoint + oSex +s(age, k=4, fx = T) + s(age, by = %s, k=4, fx = T)",
r,cv)
bv <- gamm(as.formula(bv_formula),
random = list(bblid=~1),
data = thisTable)
vc <- gamm(as.formula(vc_formula),
random = list(bblid=~1),
data = thisTable)
bic<-BIC(bv$lme,vc$lme) # get BIC
bestmod <- gsub(row.names(bic)[which.min(bic$BIC)],pattern = "$lme",replacement = "", fixed = T) #best is min BIC
#confirm there are no concurvity issues if the VC model is the best.
c<-as.data.frame(concurvity(vc$gam))
if (c["observed",contains(vars = names(c),"s(age):")]>.75) {
bestmod<-"bv"
}
switch (bestmod,
"bv" = {model <- bv},
"vc" = {model <- vc}
)
model_formula <- model$gam$formula
cat("\n\nbest model is\n",deparse(model_formula))
# Now check if the interaction is significant using parametric bootstrap `(pbootstrap=T)`
gamm_model(thisTable,
model_formula,
this_label = sprintf("%s %s final model",r,cv),
smooth_var = "age",
int_var = cv,
group_var = "bblid",
pbootstrap = T,
model_test = F)
}
}
Results for Putamen_h
NAR_Overall_Efficiency
Results for Putamen_h NAR_Overall_Efficiency M.E.
convurvity
|
|
para
|
s(age)
|
s(NAR_Overall_Efficiency)
|
|
worst
|
0.9718732
|
0.4416904
|
0.4008970
|
|
observed
|
0.9718732
|
0.3401784
|
0.1834552
|
|
estimate
|
0.9718732
|
0.3428022
|
0.3461640
|
Regression table from gamm in Putamen_h NAR_Overall_Efficiency M.E., BIC = 3412.52, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.9640954541874
|
0.19107845249858
|
88.780787327726
|
0
|
|
timepoint.L
|
0.153197738963919
|
0.495492846126982
|
0.309182544533971
|
0.757243662505699
|
|
timepoint.Q
|
0.18618333682846
|
0.36560739398768
|
0.509243904500285
|
0.610687663005984
|
|
timepoint.C
|
0.0824575819978033
|
0.178512036347206
|
0.461916090842318
|
0.644236637665198
|
|
oSex.L
|
-0.0593799783407271
|
0.056024270223721
|
-1.05989739988769
|
0.289433639425036
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
3
|
3
|
57.9559837154165
|
1.099010468818e-34
|
|
s(NAR_Overall_Efficiency)
|
2.99999999999999
|
2.99999999999999
|
2.2066659969151
|
0.0856957516751655
|
Comparing interaction models…
best model is Putamen_h ~ timepoint + oSex + s(age, k = 4, fx = T) + s(age, by = NAR_Overall_Efficiency, k = 4, fx = T)
Results for Putamen_h NAR_Overall_Efficiency final model
convurvity
|
|
para
|
s(age)
|
s(age):NAR_Overall_Efficiency
|
|
worst
|
0.9728761
|
0.5454464
|
0.5366853
|
|
observed
|
0.9728761
|
0.4972131
|
0.4513417
|
|
estimate
|
0.9728761
|
0.4943173
|
0.4777880
|
Regression table from gamm in Putamen_h NAR_Overall_Efficiency final model, BIC = 3411.55, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.8917066037838
|
0.194799823394821
|
86.7131515286214
|
0
|
|
timepoint.L
|
0.0818205246651944
|
0.495215414576723
|
0.165222087715361
|
0.868800833489742
|
|
timepoint.Q
|
0.126614359711642
|
0.365095370827026
|
0.346798041905684
|
0.728812208092882
|
|
timepoint.C
|
0.0596995807305022
|
0.178127691984599
|
0.335150475848886
|
0.737578205452899
|
|
oSex.L
|
-0.057236626890966
|
0.0557542527912067
|
-1.02658764175911
|
0.304849989212817
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999998
|
2.99999999999998
|
35.5854308842152
|
4.96099444686161e-22
|
|
s(age):NAR_Overall_Efficiency
|
4.00000000000001
|
4
|
3.66102933706076
|
0.00571354546138863
|
### NAR_F1_Social_Cognition_Efficiency
Results for Putamen_h NAR_F1_Social_Cognition_Efficiency M.E.
convurvity
|
|
para
|
s(age)
|
s(NAR_F1_Social_Cognition_Efficiency)
|
|
worst
|
0.9720741
|
0.3332461
|
0.2927798
|
|
observed
|
0.9720741
|
0.2979737
|
0.1593748
|
|
estimate
|
0.9720741
|
0.2890377
|
0.2556842
|
Regression table from gamm in Putamen_h NAR_F1_Social_Cognition_Efficiency M.E., BIC = 3415.17, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.9499019148066
|
0.191607627263226
|
88.4615198095494
|
0
|
|
timepoint.L
|
0.115240583524021
|
0.496366156549519
|
0.232168494977809
|
0.816452130959866
|
|
timepoint.Q
|
0.156982429610969
|
0.365153910655281
|
0.429907567823274
|
0.667350590496855
|
|
timepoint.C
|
0.0704308386708812
|
0.178203471409482
|
0.395227085722942
|
0.692755126412335
|
|
oSex.L
|
-0.063230721964161
|
0.056222407703891
|
-1.12465339971175
|
0.260991503852891
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999999
|
2.99999999999999
|
63.6969285137938
|
1.1090151988442e-38
|
|
s(NAR_F1_Social_Cognition_Efficiency)
|
2.99999999999999
|
2.99999999999999
|
1.3180497327979
|
0.267056369543889
|
Comparing interaction models…
best model is Putamen_h ~ timepoint + oSex + s(age, k = 4, fx = T) + s(age, by = NAR_F1_Social_Cognition_Efficiency, k = 4, fx = T)
Results for Putamen_h NAR_F1_Social_Cognition_Efficiency final model
convurvity
|
|
para
|
s(age)
|
s(age):NAR_F1_Social_Cognition_Efficiency
|
|
worst
|
0.9731073
|
0.4663956
|
0.4687256
|
|
observed
|
0.9731073
|
0.4468037
|
0.4092586
|
|
estimate
|
0.9731073
|
0.4374676
|
0.3870572
|
Regression table from gamm in Putamen_h NAR_F1_Social_Cognition_Efficiency final model, BIC = 3415.50, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.8768764482188
|
0.196003708747567
|
86.1048831986876
|
0
|
|
timepoint.L
|
0.0161656448007594
|
0.498023394943542
|
0.0324596092570953
|
0.974111665897196
|
|
timepoint.Q
|
0.080606449502637
|
0.366362092161434
|
0.220018531467274
|
0.825899287244041
|
|
timepoint.C
|
0.0453156464360474
|
0.178435390823321
|
0.253961090493068
|
0.799575103760973
|
|
oSex.L
|
-0.0657923499915266
|
0.0559581708619602
|
-1.17574161160174
|
0.239963397974342
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999998
|
2.99999999999998
|
44.1798043442418
|
3.77479651083285e-27
|
|
s(age):NAR_F1_Social_Cognition_Efficiency
|
4.00000000000001
|
4
|
2.66249153748449
|
0.0313598320078323
|
### NAR_F2_Complex_Reasoning_Efficiency
Results for Putamen_h NAR_F2_Complex_Reasoning_Efficiency M.E.
convurvity
|
|
para
|
s(age)
|
s(NAR_F2_Complex_Reasoning_Efficiency)
|
|
worst
|
0.9718918
|
0.2945251
|
0.1786761
|
|
observed
|
0.9718918
|
0.2717048
|
0.1480075
|
|
estimate
|
0.9718918
|
0.2640843
|
0.1642319
|
Regression table from gamm in Putamen_h NAR_F2_Complex_Reasoning_Efficiency M.E., BIC = 3413.37, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.9665792391014
|
0.190883357298468
|
88.8845391197317
|
0
|
|
timepoint.L
|
0.150002463671327
|
0.495006199006654
|
0.30303148520633
|
0.761925551935103
|
|
timepoint.Q
|
0.15178802091587
|
0.364976616915221
|
0.415884234444336
|
0.677579216041555
|
|
timepoint.C
|
0.0517823773328458
|
0.178094360223439
|
0.290758097380957
|
0.77129342758053
|
|
oSex.L
|
-0.0469108688629576
|
0.056271822565255
|
-0.833647582829894
|
0.40466812503491
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999999
|
2.99999999999999
|
64.8065835883413
|
2.46444471819691e-39
|
|
s(NAR_F2_Complex_Reasoning_Efficiency)
|
2.99999999999999
|
2.99999999999999
|
1.91742380908008
|
0.124989603505959
|
Comparing interaction models…
best model is Putamen_h ~ timepoint + oSex + s(age, k = 4, fx = T) + s(age, by = NAR_F2_Complex_Reasoning_Efficiency, k = 4, fx = T)
Results for Putamen_h NAR_F2_Complex_Reasoning_Efficiency final model
convurvity
|
|
para
|
s(age)
|
s(age):NAR_F2_Complex_Reasoning_Efficiency
|
|
worst
|
0.9721947
|
0.3981731
|
0.3521673
|
|
observed
|
0.9721947
|
0.3528812
|
0.1629076
|
|
estimate
|
0.9721947
|
0.3502766
|
0.2384207
|
Regression table from gamm in Putamen_h NAR_F2_Complex_Reasoning_Efficiency final model, BIC = 3410.41, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.9686246388585
|
0.191313283298937
|
88.6954859916555
|
0
|
|
timepoint.L
|
0.208327284482123
|
0.492925661486224
|
0.422634284962958
|
0.672648231259325
|
|
timepoint.Q
|
0.204636000743215
|
0.36413299562187
|
0.561981482600156
|
0.574248006732008
|
|
timepoint.C
|
0.0815260287408567
|
0.177770836781702
|
0.45860181690526
|
0.646614587991182
|
|
oSex.L
|
-0.0477755852283671
|
0.0560424609832607
|
-0.852489066149989
|
0.394136187862349
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999999
|
2.99999999999999
|
53.4014304671236
|
1.30958815159932e-32
|
|
s(age):NAR_F2_Complex_Reasoning_Efficiency
|
4.00000000000002
|
4
|
3.93386711389326
|
0.00355165292088365
|
### NAR_F3_Executive_Efficiency
Results for Putamen_h NAR_F3_Executive_Efficiency M.E.
convurvity
|
|
para
|
s(age)
|
s(NAR_F3_Executive_Efficiency)
|
|
worst
|
0.971778
|
0.5079742
|
0.4507879
|
|
observed
|
0.971778
|
0.4035207
|
0.0441854
|
|
estimate
|
0.971778
|
0.3980152
|
0.3953887
|
Regression table from gamm in Putamen_h NAR_F3_Executive_Efficiency M.E., BIC = 3414.82, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.9993746580928
|
0.190890175000003
|
89.0531671317947
|
0
|
|
timepoint.L
|
0.23111547506227
|
0.495062818286763
|
0.466840704907064
|
0.640710150331751
|
|
timepoint.Q
|
0.176573583579899
|
0.365000365354903
|
0.483762758451515
|
0.62865454599205
|
|
timepoint.C
|
0.0809055194724956
|
0.178176125023125
|
0.454076097243641
|
0.649867365965702
|
|
oSex.L
|
-0.0579282016099435
|
0.0561123005394339
|
-1.03236190733676
|
0.302138981896942
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999998
|
2.99999999999998
|
56.5263711013206
|
1.86581701703375e-34
|
|
s(NAR_F3_Executive_Efficiency)
|
3
|
3
|
1.43328545448684
|
0.231489523013777
|
Comparing interaction models…
best model is Putamen_h ~ timepoint + oSex + s(age, k = 4, fx = T) + s(age, by = NAR_F3_Executive_Efficiency, k = 4, fx = T)
Results for Putamen_h NAR_F3_Executive_Efficiency final model
convurvity
|
|
para
|
s(age)
|
s(age):NAR_F3_Executive_Efficiency
|
|
worst
|
0.9724708
|
0.5988118
|
0.5944483
|
|
observed
|
0.9724708
|
0.5273855
|
0.3925526
|
|
estimate
|
0.9724708
|
0.5190668
|
0.5110341
|
Regression table from gamm in Putamen_h NAR_F3_Executive_Efficiency final model, BIC = 3415.52, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.9290632726131
|
0.19321701809188
|
87.6168333400262
|
0
|
|
timepoint.L
|
0.140778253634304
|
0.495295678771188
|
0.284230732607177
|
0.776289340268548
|
|
timepoint.Q
|
0.124833744857213
|
0.365042196094913
|
0.341970726104102
|
0.73244103063904
|
|
timepoint.C
|
0.0463001118016715
|
0.178173654529744
|
0.259859472063207
|
0.795022876437736
|
|
oSex.L
|
-0.0484352577613098
|
0.0559378143223821
|
-0.865876837485402
|
0.386754568021397
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
3.00000000000002
|
3
|
35.4710891391769
|
7.65669845853247e-22
|
|
s(age):NAR_F3_Executive_Efficiency
|
4.00000000000003
|
4
|
2.65369901398951
|
0.0318228783401142
|
### NAR_F4_Memory_Efficiency
Results for Putamen_h NAR_F4_Memory_Efficiency M.E.
convurvity
|
|
para
|
s(age)
|
s(NAR_F4_Memory_Efficiency)
|
|
worst
|
0.9717821
|
0.3104867
|
0.1632110
|
|
observed
|
0.9717821
|
0.2893542
|
0.0942233
|
|
estimate
|
0.9717821
|
0.2749873
|
0.1272475
|
Regression table from gamm in Putamen_h NAR_F4_Memory_Efficiency M.E., BIC = 3416.06, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.979921614825
|
0.190600610952492
|
89.0863965753883
|
0
|
|
timepoint.L
|
0.189595426522223
|
0.494290502052286
|
0.383570846971621
|
0.70137375975867
|
|
timepoint.Q
|
0.185295269664508
|
0.365232732315672
|
0.507334784836196
|
0.612025867086226
|
|
timepoint.C
|
0.0848171144290412
|
0.178593463950647
|
0.47491723690672
|
0.634944168066508
|
|
oSex.L
|
-0.0564737547009858
|
0.0561878790386208
|
-1.00508785288316
|
0.315084795721971
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
3
|
3
|
65.2396514315123
|
1.05751050405678e-38
|
|
s(NAR_F4_Memory_Efficiency)
|
2.99999999999999
|
2.99999999999999
|
1.02002372858875
|
0.382902265176225
|
Comparing interaction models…
best model is Putamen_h ~ timepoint + oSex + s(age, k = 4, fx = T) + s(age, by = NAR_F4_Memory_Efficiency, k = 4, fx = T)
Results for Putamen_h NAR_F4_Memory_Efficiency final model
convurvity
|
|
para
|
s(age)
|
s(age):NAR_F4_Memory_Efficiency
|
|
worst
|
0.9721994
|
0.4696706
|
0.4041636
|
|
observed
|
0.9721994
|
0.3688750
|
0.2577590
|
|
estimate
|
0.9721994
|
0.3606478
|
0.2595430
|
Regression table from gamm in Putamen_h NAR_F4_Memory_Efficiency final model, BIC = 3422.00, obs = 1067
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
16.9692774138662
|
0.190926870512715
|
88.87841385708
|
0
|
|
timepoint.L
|
0.157915720377196
|
0.492180445582302
|
0.320849236889866
|
0.748388188357601
|
|
timepoint.Q
|
0.137889646476784
|
0.363224592722462
|
0.379626405368825
|
0.70429915851729
|
|
timepoint.C
|
0.0580375607786316
|
0.177290740625206
|
0.327358104399392
|
0.743461993685597
|
|
oSex.L
|
-0.054846163255636
|
0.0562520255977639
|
-0.975007791680592
|
0.329779897724949
|
|
term
|
edf
|
ref.df
|
statistic
|
p.value
|
|
s(age)
|
2.99999999999999
|
2.99999999999999
|
57.4095808392361
|
5.64110049193901e-35
|
|
s(age):NAR_F4_Memory_Efficiency
|
4.00000000000001
|
4
|
1.02767376286287
|
0.391772212974838
|

cog_vars <- c("NAR_Overall_Efficiency")
for (r in rois) {
cat(sprintf("\n## Results for %s\n",r))
for (cv in cog_vars) {
cat(sprintf("\n### %s\n",cv))
thisTable <- dataTable
## The logic here is to compare a bivariate interaction that varies by sex to a nested model that does not vary by sex
cat('\n### Comparing interaction models...\n')
# Bivariate interaction that just includes sex as a covariate
bv1_formula <- sprintf(
"%s ~ timepoint + oSex + te(age,%s, k=4, fx = T)",
r,cv,cv,cv)
# Bivariate that varies as a function of sex
bv2_formula <- sprintf(
"%s ~ timepoint + oSex + te(age,%s, k=4, fx = T) + te(age,%s, by = oSex, k=4, fx = T)",
r,cv,cv)
bv1 <- gamm(as.formula(bv1_formula),
random = list(bblid=~1),
data = thisTable)
bv2 <- gamm(as.formula(bv2_formula),
random = list(bblid=~1),
data = thisTable)
bic<-BIC(bv1$lme,bv2$lme) # get BIC
bestmod <- gsub(row.names(bic)[which.min(bic$BIC)],pattern = "$lme",replacement = "", fixed = T) #best is min BIC
bestmod<-"bv2"
switch (bestmod,
"bv1" = {model <- bv1},
"bv2" = {model <- bv2}
)
model_formula <- model$gam$formula
cat(sprintf("\n\nbest model is %s\n",deparse(model_formula)))
a <- anova(model$gam)
print(kable(a$pTerms.table)%>%kable_styling(position = "left"))
print(kable(a$s.table)%>%kable_styling(position = "left"))
# Now check if the interaction term from the chosen model is significant
gamm_model(thisTable,
model_formula,
this_label = sprintf("%s %s final model",r,cv),
smooth_var = "age",
int_var = cv,
group_var = "bblid",
pbootstrap = F,
model_test = F)
}
}